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SUMMARY 


Several numerical difficulties arise in the problem of integrating the 
ordinary differential equations that represent the flow of a gas in chemical 
nonequilibrium through a nozzle. These difficulties are identified with par- 
asitic eigenvalues, saddle-point behavior, and indeterminate forms. The 
numerical difficulties and the methods that successfully overcame them for the 
particular examples analyzed are discussed in this report. 


INTRODUCTION 


The equations representing one -dimensional steady flow of a gas in chem- 
ical equilibrium through a nozzle (channel flow) are well known (see, e.g., 
ref. l) . These equations are treated again here for the purpose of studying 
certain difficulties, specifically, saddle-point behavior, that arise when 
they are integrated numerically through the transonic region. The technique 
developed from the simpler problem is then applied to clarify some of the dif- 
ficulties encountered in integrating the much more complicated equations that 
represent the flow through a nozzle of a gas in nonequilibrium. The physical 
model of air used in these computations is the same as that used in refer- 
ence 2. This model has been greatly improved in the ensuing years in a vari- 
ety of ways by a variety of authors. However, this report is an analysis of 
numerical difficulties that can arise in integrating the equations formed 
from the choice of any model, and, for this purpose, the one chosen appears 
reasonably representative. 

In reference 2, a critical analysis of several numerical methods for 
computing the flow of a gas in chemical nonequilibrium has been carried out. 
Certain information contained in that reference, such as the concept of para- 
sitic eigenvalues and the construction and value of implicit methods, is help- 
ful for a detailed understanding of the material presented here. However, in 
reference 2, the examples are limited to the nonequilibrium flow behind a 
normal shock wave. Additional complications arise in the transonic section of 
a nozzle; two of the most important are those arising from the existence of a 
saddle point and the problem of numerical solution of the differential equa- 
tions when the gas they represent is nearly in equilibrium. 



The question of integrating through the saddle point has generally been 
avoided by considering the gas to be in equilibrium from a position upstream 
of the minimum section to one a little downstream from it. In fact, most 
often a pressure distribution along the channel is assumed in this region; 
the corresponding cross-sectional area follows from a direct calculation, and 
a few iterations produce the desired approximate nozzle shape in the transonic 
part. This procedure is not always usable, however, and, in this report, 
direct methods are developed for treating saddle points. For example, if the 
flow is out of equilibrium upstream of the throat, the abovementioned device 
is unsatisfactory. An example of such a case is included here. Also, in 
other problems, such as blunt -body flow, the method of integral relations 
often leads to saddle-point behavior, and it becomes necessary to use 
techniques such as those to be presented. 


SYMBOLS 


[ ] 

[ r 

A 

Ax 

[A n ] 



matrix of enclosed quantity 

inverse of matrix 

cross-sectional area of nozzle 

derivative of A with respect to x 

matrix in locally linearized equations (see eq. (10)) 

equilibrium speed of sound 

critical speed of sound 

element in A n matrix 

?n c P ± 

specific heat at constant pressure for ith species 


er A 

F 

h 


[I] 


n 


truncation error 

derivative of w with respect to s (see eq. (6)) 
step size 

enthalpy of ith species 


unit matrix 
step number 
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p pressure 

production of species i in moles per unit volume per unit time 
B universal gas constant 

s independent variable 

s Q scale factor (see eqs . (3)) 

T temperature 

u velocity 

w dependent variable in coupled equations (see eq. (6)) 

x distance along nozzle 



7 ratio of specific heats for a perfect gas 

7^ molar concentration of ith species, moles/gm 

p density 

a eigenvalue of matrix [A n ] 

Superscripts 

vector 

T transpose of vector 

* differentiation with respect to s 

BASIC EQUATIONS 

Equilibrium Flow 

When an inviscid channel flow is in chemical equilibrium, the basic 
equations can be expressed in terms of the three dependent variables u, p, 
and p by the matrix equation 
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where x is the independent variable representing length along the channel, A 
is the channel cross-sectional area at a given x, and the subscript indicates 
differentiation with respect to x. It is well known that the solution of 
equation (l) can be found directly^ and the dependent variables u, p, and p 
can be expressed as functions of A and reference conditions (see, e.g., 
ref. l) . As was pointed out in the Introduction, however, we are interested 
here in the numerical solution of equation (l) in order to study the tech- 
niques required to treat more complicated problems with fundamentally similar 
numerical difficulties that cannot be resolved analytically . 


Equation (l) can be solved directly for the derivative terms. One finds 

_ \ 

ua 2 A x 

^ A(a 2 - u 2 ) 


p x n 


A (a 2 


pu 2 A x 
2 - u 2 ) 


p^a 2 ^. 

Px = A (a 2 - u 2 ) 

J 

which give, by inspection, the well-known condition that if the velocity 
becomes sonic, the derivatives of the dependent variables all have zero denom- 
inators; if they are to remain continuous, the numerators must also vanish. 

In this case, the numerators can vanish only if A x = 0, which leads to the 
well-known fact that the sonic velocity occurs at the throat . Of principal 
interest for our purposes is the fact that this indeterminate form leads to 
the study of critical points (in particular, saddle points) and the numerical 
difficulties that occur near them. 


An analysis of saddle points from a numerical standpoint is presented in 
the next part. Now we note only that the indeterminacy in equations ( 2 ) can 
be eliminated by a sinple transformation; the general procedure is outlined 


in reference 2. Introducing the new independent variable 


defined by 


x* = = SoA.(a 2 - u 2 ) 

where s Q is a constant scale factor, makes x a dependent variable in the 
new set of four coupled nonlinear differential equations 

u' = ff = -s 0 ua 2 A x (3a 

p' = s 0 pu 2 A x (3t 

p' = s Q pu 2 a 2 A x (3c 


x* = s Q A(a 2 - u 2 ) 


in ■ vi 1 in ■■ il 1 ml 11 1 1 


Ilf 



Notice that equations (3) are determinate for all finite values of s and 
autonomous (i.e., s does not appear explicitly on the right-hand sides of the 
equations) . 


Nonequilibrium Flow 


The basic equations used herein for nonequilibrium channel flow are 
identical to those presented as equations (l) in reference 2. With a slight 
change in the notation (z* for ¥>y±> Cp ^ or ?^i°p-j_^ ’ equations can be 


written 




t- du t- d.p 
pA. — + uA. — 
^ dx dx 


pu 


du 

dx 


N 

„ dp * dT V d ^i 

+ RzT a£ + ^ z sr 


= -pu/U 


= 0 


u 


du 

dx 


_* dT 

+ Go — + 
P dx 


N 


, d 7i n 

h,* 3 — = 0 
1 dx 


pu 


d7x 

dx 


= Q 1 (p,T,7i, . . 7 n ) 


(4) 


N N / \ 

pa — = Q ?y ± > • * *9 7^) 

If equations (4) are solved for the derivatives du/ dx, dp/ dx, and dl/dx, 
each expression has in its denominator the term 

TRz*Cp - u 2 (C* - Rz*) 

which is the counterpart of the term a 2 - u 2 that appeared in equations ( 2 ) . 
It can be eliminated by a transformation similar to equation (3&) • There 
results 
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s 0 puRA^-z*Di ^ 

CpTD 2 

- Rz*T P u 2 cJa x 


s 0 p a RA’|z*D 1 - 

c p td J 

+ p 2 u^Cp - Rz*) 

Ax 

s oP a[(u 2 - Rz* 

r T)D 1 - 

u 2 pRTD 2 ] + Rz*Tpu s A x 

s 0 puA[TRz*Cp ■ 

■ u2 C c p 

- Bz *)] 


s oP-A-Q ' [j® z Cp 

- u 2 (c* 

- Rz *)] > 1 = 

1, 


where 






N 


J 


and where s D is again a scale constant. The dependent variable T has 
replaced p, used in the equilibrium case, and the N species concentrations 
are now coupled into the equations. Notice, however, that the equations are 
determinate and autonomous. 

Some of the numerical difficulties that occur in integrating equa- 
tions (5) can be illustrated by a study of the much simpler equations ( 3 ). 
Hence, the first part of the report will be devoted to the numerical analysis 
of equilibrium flow. 


Local Linearization 


The locally linearized forms of equations (3) and (5) will next be con- 
structed, following the procedure of reference 2. In either case, equilibrium 
or nonequilibrium, the equations of motion can be represented in vector 
notation as 


where, for equations ( 3 ), 


_ dw _ 

w = — = F(,w; 
ds 


v 1 = (u,p,p,x) 


(6) 

(7) 
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and for equations (5) 


• 9 7;pp x ) 


( 8 ) 




- (u,p ,*£, 7^9 


Note that equations (6) are quas i -linear , that is, the highest derivative 
appears linearly. The linearization is now performed by expansion of each 
component of the vector F, denoted Fp, in a multidimensional Taylor series 
about the point s n = nh. This gives, a:ffcer substituting equation (6), 



for the ith component w^_ of w f . The number of components of w is taken 
to be m; hence, there are m coupled equations such as ( 9 ) in the set. Note 
that when w - w(s + nh) = w n+1 


I w - w. 


nl 2 = Rn+i “ ^nl 2 53 ^1 ^ + °( h3 ) 


by Taylor expansion. By defining the Jacobian matrix 


[A] = (a id ) 



the set of equations (6) can be written 

w’ = F n + [A n ](w - w n ) + 0(h 2 ) 


= [A n ]w + f n + 0(h 2 ) (lO) 

where 

?n “ F n - [A n ]w n (11) 

It is important that although equation (lO) represents the local linear- 
ization of equation (6) (referenced to s n = nh) , it is the derivative w 1 
that is represented as linear in the stepsize h. Therefore, in a differenc- 
ing scheme, the linearization here is consistent with any method that approx- 
imates the function itself with a second-degree polynomial. For example, if 
the difference scheme is 


w n+i = w n + 2 


h ^n+i + *%) + 0(h 3 ) 
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then inserting expression (lO) for the derivative terms does not alter the 
order of the error . 

The [A] matrix is not difficult to construct analytically in the case of 
equilibrium flow represented by equation (3) • It is especially simple (and 
just as instructive) to asstime a calorically perfect gas; then 


7P 3a 2 _ a^ 3a 2 _ aP_ 

P ’ 3p P ' 3p “ P 


and 


[A] = 


-s 0 a 2 Ax 

2s 0 puA x 

2s c pua 2 A x 

-2s 0 uA 


1 27 

So - ua A x 

s 0 u 2 A x 

0 

-s 0 ^ a 2 A 


-s 0 p ua 2 A x 


0 


“ s 0 ua 2 A xx 

s o pu ^XX 


s Q - pu 2 a 2 A x SoPU^A^ 

s 0 — a 2 A s (a 2 - u 2 )A 

u p o' ' x 


( 12 ) 


If the nonequilibrium equations (4) are to be integrated by means of 
implicit formulas, the A matrix will be 13 x 13 for the gas model used in 
reference 2 and in the present report. However, just as in reference 2, the 
elements of A are not determined analytically. They are computed by the 
approximate formula 


3P t F i (l.01w j ) - F 1 (0.99wj) 

3wj 0.02w^ 


(13) 


NUMERICAL INTEGRATION THROUGH A SADDLE POINT 


Introductory Example 

The analytical behavior of autonomous systems of first-order ordinary 
differential equations has been thoroughly studied in the mathematical litera- 
ture. The facet of this material of particular interest in the present inves- 
tigation is well summarized in reference Z, page 70, ff. As mentioned 
previously, this part of the theory is the study of behavior of solutions near 
a critical point, specifically, a critical point of the type known as a saddle 
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point. Because of the advantage of geometrical visualization of the points 
involved, the case involving only two independent variables will be considered 
first . 


Let there be two differential equations 


if * p ( x ’”> 1 
f ■ ^ 


which can be put in the form 


c 3x 


= aai x + a 12 w + f x (x,w) 


dw 


uw ■ , „ / \ 

— = a 21 x + a 22 w + f 2 (x,w) 

Further , suppose that a 11 a 22 - a 21 a 22 / 0, and that 

f n 


lim 


x,w -»0 */x 2 + w 2 


= 0 , i = 1, 2 


( 14 ) 


(14a) 


(l4b) 


The point x = w = 0 is then said to be a simple singularity of the system. 
It is shown in the theory (see, e.g., ref. 3) that under such circumstances, 
the singularity of system (l4) at x =• w = 0 has essentially the same nature 
as that of the linearized system 


dx 

dt 


a 1:L x + a 12 w 


dw 

dt 


3<pi x “f - a P p) W 


(15) 


The constants a** are components of the Jacobian matrix used above, that is. 


a 


li 




Hence the condition on the constants a-^ can be stated in terms of the non- 
vanishing of the Jacobian S(P,Q)/S(x, w) at the origin. 

Note that in this case the differential equations (l5) can be combined to 

give 

dx = a xl x + a lg w 
dw a 21 x + a 22 w 
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The point x =• w =■ 0 is then seen to be one at which the derivative is indeter- 
minate. While this simple means serves to fix the saddle point when there are 
two variables involved, such as x and w above, it is not clear at once how 
it can be extended to more complicated cases, such as that represented by the 
system of equations ( 5 )* A more suitable test for more than two variables is 
given below. 


As an example, consider the autonomous linear system of differential 
equations : 


or 



x + 0-5 



+ 1.0 


dw _ 2w - x + 0.5 
dx 1 - w 


( 16 ) 


(l6a) 


From the latter form it is seen that the critical point is at w = 1, x = 2-5- 
It is helpful to put equations (l 6) into the vector-matrix notation used 
above in the section on Local Linearization. Thus, put 


W 1 = (w,x); 


[A] 


~ 2.0 - 1 . 0 " 
= 1.0 0 J ; 


f T = (0.5, 1.0) 


and equations (l6) can be written 


If’ = [A] w + f 

This equation is formally equivalent to the linearized form, that 
obtained by neglecting terms 0(h 2 ) in equation (10) . 

The eigenvalues of the matrix in equation (l6b) are 

a x =1 - J2 « -O.blk 

cr 2 = 1 + s/2 « 2.in4 


(l6b) 


is, that 


( 17 ) 


They are real and of opposite sign. This fact is important since it character- 
izes, for the two-dimensional case, the critical point as a saddle point. In 
cases where more variables are present, the presence of eigenvalues of differ- 
ent sign, in company with a requirement that solutions remain finite, can be 
taken as indications of saddle-type behavior. The behavior of the integral 
curves of equation (l6a) (or (l6b)) is easily found. Independent solutions 
for w and x are combinations of the complementary solutions e cTlS and e 02 
with the particular integral w = 1.0, x = 2-5- If and x Q are the values 
of w and x, respectively, at s =• 0, then the general solution of equa- 
tion (l6b) is 

10 


I 



w - 1.0 


1 


2 \]~2 


[ (x D - 2.5) - <y ± (w Q - 1) ] 


ffis >. 


^j= [(x c - 2-5) - cr 2 (wo - l) ]e° 2S 


x - 2*5 ~ 


^2 

2j2 


[(x Q - 2-5) - CTi(w 0 - 1 )] 


01 S 


2 s/2 


[(x 0 - 2 . 5 ) - o 2 (w 0 - l)]e GsS J 


(18) 


Since and cr 2 are real and of opposite sign, with cr 2 > 0, it is 

clear that one or the other of the exponential terms must have a zero coeffi- 
cient if a solution is to remain finite as s -> ±oo. Thus, a finite solution 
as s -> +oo results if x G - 2-5 = cr 2 (w 0 - l) while a finite solution as 
s -» «oo results if x G - 2. 5 = cri( w q - l) . In these two cases, the integral 
curves are contained in the lines (separatrices) 


w 



2-5) 


( 19 a) 


W 



2-5) 


( 19 b) 



Sketch (a).- Saddle curve behavior for 
equations (l6). 


respectively. These lines are 
identified in sketch (a) . 

If the initial conditions (x Q , 
w G ) lie on line (2) in the sketch, 
then as s increases the solution 
will follow this line toward the 
point (2-5* 1.0). Integration in 
the opposite direction (s decreas- 
ing) would give an integral curve 
moving away from the saddle point. 
Similar remarks, with sign reversal 
on s, hold for the line (l) in the 
sketch. 
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If the initial conditions do not satisfy one of equations (19), no 
matter how small the error , the corresponding integral curve will contain an 
exponentially increasing component that will prevent the curve from ever 
approaching the saddle point. Instead, the curve will have a hyperbola -like 
shape and will ultimately approach ±°°. This fact is at the root of numerical 
difficulties in the solution of differential equations having saddle-point 
singularities; since computing machines do not reproduce numbers exactly, 
there will always be an error in the initial conditions. 


We shall now discuss the numerical determination of the integral curves 
of equation (l6b). The particular method chosen is the second -order Runge- 
Kutta: 


(i) 

<+i = % + hu n 

u n+i = u n + I h ( u n+l + \) 


( 20 ) 


This method will also be used later in integrating the flow equations dis- 
cussed in the previous sections. Precise meanings of terms and symbols used 
in equations (20) and the numerical analysis that follows can be found in 
references 4- and 5- 


Accuracy 

The accuracy of the integral curves found by applying the integration 
method (20) to the differential equations (l 6) will be discussed first. The 
truncation error, er^, for the second -order Runge-Kutta method is given by 
(see ref. 5, p- 63) 

er A = \ ( CTh ) 3 


Here, a = cr 2 = 2.4l is the largest eigenvalue. Hence, a stepsize 
h = As = 0.05 gives a local truncation error of about 

er A = 3x10 " 4 

Thus, we should expect a total accumulated ("global") error of less than 
3 percent after 100 steps. This expectation is borne out by comparison of 
the calculated solution with the exact one. The results (which correspond to 
curves AA* and AA' 1 in sketch (a)) are listed below (note that the initial 
conditions were chosen to lie above and below the separatrix by approximately 
equal amounts) . 
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TABLE I 


Curve 

Initial values 

Numerical results 
after 100 steps 

Exact results 
for same s 


w o 

x o ■ 

w 

X 

w 

X 

AA’ 

-0.035531 

0 

1.2411 

2.0232 

1.2516 

2.0189 

AA’ » 

-0.035537 

0 

.4686 

2.3413 

.4595 

2.3470 


The point to be made here is that the numerical method ( 20 ) is suffi- 
ciently accurate, that is, increasing the accuracy of numerical integration 
would not lead to discernibly different results to the scale of sketch (a) . 


Instabilities 

Induced instability . - The real stability boundary of the Runge-Kutta 
method given by equation ( 20 ) is - 2.0 (ref. 5 , p. 83). That is, in the 
present case, h is limited to the range 

| ah | <2 

where a is the largest negative eigenvalue of the matrix [A] in equa- 
tion (l 6 b). This value is a ± = - 0 - 4 l 4 , so that h can be as large as 4.8 
without bringing on the induced instability. If the direction of integration 
were reversed, the largest negative eigenvalue would become -cr 2 - - 2 . 4 l, and 
there would be no induced instability for h < O.83. Hence, for the present 
problem, there is absolutely no danger of induced instability with the 
stepsize h = 0.05 used. 


In herent instability . - Prom the above discussion, it is seen that there 
is actually no numerical difficulty, either in regard to accuracy or to 
induced instability, in calculating the integral curves marked AA T and AA ! T 
in sketch (a), the initial values for which are given in table I. Thus, if 
such a result is indeed the correct answer for the problem set, there is no 
more to be said. However, it is usually the case that when a saddle point 
occurs in a problem the desired result for the integral curve is the separa- 
trix or saddle curve. As we have seen from the example above (sketch (a)), 
reaching the saddle point by forward integration from a given initial posi- 
tion is difficult indeed. Por the case examined above, the initial condi- 
tions for a point on the saddle line starting at x = 0 are 


Xq = 0, W Q = | (7 - 5 s/ 2 ) 

Since the computer is unable to represent exactly the number n/~2", it is 

impossible to exclude completely terms in e°* 2S = g ( 1 ' + n/2)s ^ om solution 
(s increasing) . 
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This is an example of an inherent instability (ref. 4, p. 48). That is, 
it is not possible to exclude from the numerical solution the terms in equa- 
tions (l8) that contain the growing exponential, e 2 - e 2-4i4S a Again, this 

occurs because an irrational number cannot be represented exactly on a binary 
digital computer. 

In summary 

1. There is no numerical difficulty per se, either with regard to 
accuracy or stability, in numerically evaluating integral curves in the 
Vicinity" of a saddle curve. (The precise definition of vicinity is gov- 
erned by the word length of the particular computer used for the computation.) 

2. The saddle curve itself is inherently unstable and its numerical 
calculation involves all the difficulties associated with such an instability. 


Numerical Integration Through Saddle Points 


General dicussion .- Actually, the numerical integration of saddle curves 
may or may not lead to difficulty. Recall that, by definition, a saddle 
point can occur only if there are two real eigenvalues of opposite sign in 
the [A] matrix of the coupled equations. No matter in what direction the 
integration proceeds, therefore, there is always a positive eigenvalue in the 
presence of a negative one; and, again by definition, an inherent instability 
exists if that particular solution is desired for which all effects of the 
positive eigenvalue disappear from the exact analysis. We now show that the 
numerical difficulties depend upon the ratio of these two eigenvalues and the 
direction of integration when their magnitudes are unequal. 


Consider equations (l5) • Transform, by an appropriate rotation and dis- 
placement, to a new set of variables, y x and y 2 , that lie along the saddle 
curves, with origin at the saddle point itself (as shown in sketch (b)). 

This transformation uncouples equa- 
tions (15) and we find solutions in 
the form 



7i(s) = yi( 0 )e aiS 
y P (s) = y P (o)e^ s 


( 21 ) 


The parameter s can be eliminated 
from these equations, leading to the 
solution 


Sketch (b) 

-y x ( s)- 

l/Cfl 

y 2 (s) 


_y 1 (°)_ 


_y 2 (o)_ 


i/°l 


( 22 ) 


Now if y 1 is started a distance -2-5 units from the origin (saddle point) 
and y 2 a distance e above the axis, we have 
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Yi(s) = -2.5 


e y 2 ( s ) 


■Fi/Oa 


If we proceed to integrate (accurately) to where y 2 = 1, there results 


yi 



+^S) 


-2. 


- 0.171 

5e 


Even if e had been chosen as small as 10 " 6 , y ± would still be -0.223 units 
from the y 2 axis after the integration. This accounts for the divergence 
from the separatrix of curves AA 1 and AA* 1 in sketch (a). 

Ne:xfc, notice that if conditions are reversed and one starts with 
(y x = e, y 2 = -2.5) and proceeds to integrate until y x = 1, the results 
would be 


y? 


-2.5e 


5 .83 


This leads to a behavior quite different from that encountered in integrating 
from A to A T . To illustrate this, integration of equations (l6) was started 
at point B in sketch (a), and integrated (using eqs . (20) with As = -0.05) 
for 140 steps. The results are given by the solid line BB 1 in the sketch, 
and the initial and final values are presented in table II. 

TABLE II 


Curve Initial values 


Numerical results 
after l40 steps 


Exact results 
for same s 


BB ! 


w x 

0 2.9 


W X 

0.91060 2.28417 


W X 

0.91059 2.28414 


This difference in divergence of the numerical solution from the saddle 
point can be utilized. Suppose that one has an initial value problem in 
which the initial data are given on the saddle curve itself at a point such as 
A in sketch (a); the problem is to integrate numerically across the saddle 
point and proceed along the same saddle curve. Suppose, further, that pro- 
ceeding forward (increasing s) from A, one found, by experiment, results 
similar to those represented by curves AA T and AA* 1 . It would appear, on the 
basis of the discussion in the previous paragraph, that one could readily 
solve the problem by iterating on results obtained by integrating backward 
from suitable guesses near the critical point (but still in the same half- 
plane as the initial point). This approach does not succeed in all cases; 
there are two conditions under which it either is not helpful or fails 
disastrously: 

1. The eigenvalues cr 1 and cr 2 are equal in magnitude so that the inte- 
gral curves are right hyperbolas and neither direction of approach to the 
saddle point is preferred. 
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2. Coupled into the equations are other eigenvalues that are relatively 
large negative numbers with respect to forward integration. These become 
large positive numbers when the direction of integration is reversed and lead 
to serious inherent instability. 


It will be found that both the above conditions operate in the problem of non- 
equilibrium channel flow when passing through the throat section. Neverthe- 
less, the saddle points that arise in the solution of the blunt-body problem 
by means of the method of integral relations do indeed approach their asymp- 
totes at different rates. Thus the device mentioned above may be put to good 
use. 


Some practical considerations.' 


equations ■fl6) 

forward integration, the saddle point becomes the point at 


If each of the quantities w and x in 
is considered as a function of s (see eqs. (l8)), then, for 

s =■ +oo. The deriv- 


atives w 1 and x 1 approach zero as s 



Sketch (c).- Saddle curve behavior in terms of 
independent variable s. 


00 , as shown for w(s) in sketch (c). 
Care should be taken that too much 
calculation time is not wasted in 
indiscriminate attempts to integrate 
too close to the saddle point in 
terms of s. 

An extrapolation across the 
saddle point to obtain initial 
values from which the solution can be 
started again on the opposite side is 
often satisfactory. The extrapola- 
tion should be based on two integral 
curves that are initially very close 
and later diverge in opposite direc- 
tions, such as curves AA* and AA" in 
sketch (a) . To offset the inaccuracy 
inherent in such an extrapolation. 


notice that integral curves have the property that the direction of motion 
along the curve with increasing s reverses as the saddle point is crossed. 
Hence, by extrapolating past the saddle point, and reversing the sign of the 
increment As, an integral curve that approaches the saddle curve, rather than 
one diverging from it, will be obtained. This convergence will not occur if 
the saddle point is not passed in the extrapolation. For example, starting at 
point C in sketch (a) and integrating 134 steps with stepsize As = h = -0.05, 
the curve CC 1 was obtained. The end results are given in table III. 


TABLE III 


Curve Initial values 

w x 

CC’ 1.2 2.6 


Numerical results Exact results 
after 13 ^4- steps at same s 


w x 

2.0156 4.9518 


W X 

2.0158 4.9523 


A point on the saddle line itself occurs at w = 4*9528, x = 2.0l60, so that, 
even with the "bad start" at C, the integration approaches (for this example) 
the desired saddle curve with good accuracy. 
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NUMERICAL CALCUIATION OF EQUILIBRIUM FLOW 


The equations for inviscid channel flow of a gas in chemical equilibrium 
have been given above (eqs. (l)). These equations were then put in autonomous 
form (eqs. (3)), and they will be treated in that form in the numerical work 
described below. In the following discussion of the ‘numerical solution of 
equations (3), the eigenvalue structure at the throat is first considered. 

Then the integration process used to approach and pass through the saddle 
point at the throat of the channel is presented. 


Conditions at the Throat 


Although all solutions for equilibrium channel flows presented in this 
report are for a real gas in chemical equilibrium (carried out by means of 
the real -gas program described in ref. 6 ), it is instructive to inspect ana- 
lytically the nature of the solution at the throat when perfect gas approxima- 
tions are valid. This is because the eigenvalues (as determined numerically) 
at the throat, for the equilibrium real -gas flow, show a similarity with 
those determined analytically for a perfect gas. 


Thus, consider the locally linearized equations for perfect -gas, channel 
flow right at the throat section where A x = 0 and u = a = a*. The A 
matrix given by equation (l2) simplifies to 


0 0 


0 


-s 0 a 


* 3 a. 


■XX 


[A] = 


0 

0 


0 

0 


0 Sopa^Axx 

0 s oP a* 4 A 


-2s 0 a*A 


-s 


o 


1 

p 


a* 2 A 



0 


(23) 


Note that the only differences between this result and the corresponding one 
for a real gas appear in the elements a 42 and a43j where the terms (Sa 2 /Sp)^ 

and (5a 2 /dp)p have been simplified in the present case of a perfect gas. 


The eigenvalues of the matrix in equation (23) are given by the roots of 
the equation 


2r; 


a - s o AA-^a* 


4 ( £. a.*2 


+ 1 


= 0 


P O / p 

After some reduction, using the relation a*^ = a = j the eigenvalues 
are found to be 
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— i[^(7 + l)AA-xx a s q 2 J 


1/2 


(24) 


*3 = *4 = ° 


The scaling factor s Q , still at our disposal, is taken to be 

1 


s ° “ a 2 
a o 


(25) 


where a Q Is the sound speed at u => 0, that is, under reservoir conditions. 
This is convenient because then 

dx a 2 - u 2 r- t- 
— = — A as A 

ds o 2 
a Q 

when u « a. In all cases considered in this report, the area variation is 
taken to be 

A = 1 + x 2 (26) 

With the above choices for s Q and A, the nonzero eigenvalues become 



- O * 4 

= ± 

J 

(7 + 1)2 

a O - 


1/2 


= + 


~\7 + 1 


1/2 


±1.83 


(27) 


when 7 = 1.4. The eigenvalues at the throat in a perfect-gas flow therefore 
are real and of opposite sign. Hence, there is a saddle point at the throat, 
and the situation is one in which backward integration near (but not across) 
the saddle point becomes useless, as pointed out previously. The same situa- 
tion is found when we deal with a real gas, where the eigenvalues are 
numerically determined. 


Calculations Upstream and Downstream of the Throat 

The channel flow equations in autonomous form, equations ( 3 ), were inte- 
grated numerically by means of the second -order Runge-Kutta procedure dis- 
cussed above (eqs. (20)). The first case considered is one in which the 
approximate reservoir conditions are 
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u x I 0 ' 5 , cm /sec 



Sketch (&).- Integral curves near the saddle 
curve for a real gas flowing in equilibrium 
through a nozzle; reservoir conditions, 
equations (28a) . 



Sketch (e).- Integral curves near the saddle 
curve for a real gas flowing in equilibrium 
through a nozzle; reservoir conditions, 
equations (28a). 



Sketch (f ) . - Integral curves near the saddle 
curve for a real gas flowing in equilibrium 
through a nozzle; reservoir conditions, 
equations (28a). 


T 0 = 10,000° k 
P Q = 4470 atm - 

Integration was begun at x = -1.0 cm, 
where the throat, for the channel 
shape given by equation ( 26 ), occurs 
at x = 0. Results for four guessed 
values of the initial velocity are 
shown in sketch (d) . The saddle- 
point behavior near x = 0 shows 
clearly. Also, sketches (e) and (f) 
show the calculations of density p 
and temperature T corresponding to 
the initial guessed values of veloc- 
ity. Many such curves can be calcu- 
lated in a few seconds, and the 
computations can easily be automated. 
Thus, in the final program, such solu- 
tions are calculated until the saddle 
point is sufficiently well located, 
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an extrapolation is made across the saddle point, and the integration con- 
tinued as far as required downstream, all in a single run. The final down- 
stream variation of the physical quantities u, p, p, and T is shown in 
sketch (g) . These results will be used later in comparisons to be made with 
calculations of nonequilibrium flow. 


In support of the claim made above that the eigenvalues of the equilib- 
rium real-gas flow calculations are analogous to those for a perfect -gas 
flow, sketch (h) shows the eigenvalues actually calculated (for reservoir con 
ditions (28a)) as functions of distance x along the channel. It is seen 
that, indeed, when x = 0, two of the eigenvalues vanish and the remaining 
two are equal in magnitude and of opposite sign. (The magnitude of the eigen 
values is 2.0 rather than 1 . 83 - ) In this sketch, only the real part of the 
complex eigenvalues is shown (all complex eigenvalues appear as conjugate 
complex pairs since the elements of the matrix [A] are real) . The maximum 
magnitude of any eigenvalue is less than 6 in this case, so a step size 
h =As - 0.02 is well within the stability boundary (|crh| < 2) for the inte- 
gration method used. At no stage of the calculations were any of the eigen- 



Sketeh (g).- Channel flow solution for a real gas 
in equilibrium downstream from the throat of a 
nozzle; reservoir conditions, T Q = 10,000° K, 
p Q s 4 ^ 1-70 atmospheres. 


values pure imaginary, for which the 
Runge-Kutta second-order method is 
unstable (see ref. 5 , p- 83)- Hence 
the explicit equations (20) were 
adequate for solution throughout the 
entire range. 



Sketch (h).- Local eigenvalues in the equations 
representing the channel flow of a real gas 
in equilibrium flowing through a nozzle; 
reservoir conditions, T 0 » 10,000° K, 

P 0 - Mj- 70 atmospheres. 


20 



The eigenvalue patterns for two other solutions, corresponding to 
assumed equilibrium channel flow with reservoir conditions 


T n = 8,000 K 


p = 100 atm 


( 28 b) 


and 


T 0 = 10,000 K 


p =0.2 atm 

o 


( 28 c) 


are shown in sketches (i) and (j). The scaling (factor s Q in eqs. (3)) was 
chosen so that the right side of equation (3d) was 1 at the beginning of the 
integration. It is interesting to note that in spite of the large difference 
in reservoir pressures (0.2 to 4470 atm), the eigenvalue curves are much 
alike in shape and magnitude in the three cases. The nonequilibrium real-gas 
flow for the same three sets of reservoir conditions is examined in the next 
section. 




Sketch (i).- Local eigenvalues in the equations 
representing the channel flow of a real gas 
in equilibrium flowing through a nozzle; 
reservoir conditions, T 0 = 8000° K, 
p Q = 100 atmospheres. 


Sketch (j).- Local eigenvalues in the equations 
representing the channel flow of a real gas 
in equilibrium flowing through a nozzle; 
reservoir conditions, T 0 = 10,000° K, 
p Q = 0.2 atmospheres. 
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NUMERICAL CALCULATION OF NONEQUILIBRIUM FLOW 


Introduction 

A variety of problems can arise when one seeks to make numerical calcula- 
tions of a gas flowing out of chemical equilibrium. Some of these are dis- 
cussed in the following. The problem of providing a suitable mathematical 
model Is not considered here; we use the same model as that described in ref- 
erence 2, which, in turn, was taken from references 7 and 8. There is no 
essential difficulty in casting the equations in finite-difference form. In 
fact, there is little numerical difficulty in constructing the [A n ] matrix 
(see eq. (l3)) if it is required. The real difficulties are more fundamental. 
For example, the saddle-point problem still occurs in the throat region, and, 
if treated in the way outlined above, requires an extrapolation of all nine 
species across the critical point (which is no longer exactly at the minimum 
section) . In addition, the problem of parasitic eigenvalues (numerically 
large negative eigenvalues which, for numerical stability, force an unnatu- 
rally small step size, see ref. 2 ) can now arise. Furthermore, these para- 
sitic eigenvalues appear near the beginning of the calculations (i.e., in the 
throat section) rather than downstream as is the case in the study of 
nonequilibrium flow behind shocks. 

The above difficulties can be compounded in some calculations .when the 
gas is very nearly in equilibrium, because, in such a case, the Q 1 in equa- 
tions (4) and (5) are composed of terms formed by the product of very small 
and very large terms (indeterminate forms in the mathematical model when the 
equilibrium limit is approached) . Surprisingly, this near indeterminacy 
appeared to cause no trouble in any of the calculations carried out. This 
matter is reserved for discussion in the last section after the results have 
been presented. 


Parasitic Eigenvalues 

It is well known that the equations representing a gas flowing in chemi- 
cal equilibrium have local parasitic eigenvalues. That is, some of the eigen- 
values of the [A n ] matrix (found by putting eqs. (5) in the form of eq. (10)) 
are very large negative numbers compared to others that determine the solu- 
tion. 1 In flows from nozzles with high reservoir pressures, these parasitic 
eigenvalues appear upstream of the throat. 

Some of the eigenvalues for flow from the reservoir conditions given in 
equation (28a) are shown 2 in figure 1. Actually, several large negative 
eigenvalues are contained in the local [A n ] matrix. The three largest ones 
are shown in the figure. Now it should be mentioned that, for these reservoir 
conditions, the gas in the nozzle between the reservoir and a section about 


1 These eigenvalues, which determine the solution, are called driving 
eigenvalues . 

2 Notice the numerical evidence of singular ' behavior at the critical 
point . 


22 



3 cm downstream from the throat is very nearly in equilibrium. This can he 
shown by comparing solutions for a nonequilibrium case with those for an equi- 
librium one (fig. 2) . In this region of near equilibrium, the eigenvalues 
that are actually driving the solution are those already shown in sketch (h) , 
two of which are also included in figure 1 for comparison. This extreme dif- 
ference (about lxlO r ) between the largest negative eigenvalues and the driving 
eigenvalues forces one, practically speaking, to an implicit numerical 
integration method. 

As in reference 2 , equations ( 5 ) were converted to equations (10) by cal- 
culating the elements ( a ij) n numerically, that is, by using equation ( 13 ). 

The modified Euler implicit differencing scheme 

w n+i = w n + |h(v^ +1 + v£) ( 30 ) 

was used, leading to the formula 

([I] - | h[A n ]) (w n+1 - w n ) = hF n + 0(h 3 ) (31) 

Calculations were started at x = -1 cm, and all integrations were carried 
out with s as the independent variable. Although the product of the step 
size and the largest negative eigenvalue started as high as - 8000 , the results 
were satisfactory, as illustrated in figure 2, where the nonequilibrium cal- 
culations are compared with the previously discussed equilibrium ones. The 
equilibrium species concentrations shown in figure 2(b) are readily calcu- 
lated for the given reservoir composition from the assumption of chemical 
equilibrium and knowledge of two thermodynamic variables, say, density p 
and temperature T as shown in figure 2(a). 

The saddle -point problem was treated exactly as in the equilibrium case 
by finding two integral curves that diverged to opposite sides of the criti- 
cal point. The results for p, u, and T are shown in sketches (k), (& ) , and 
(m) . Extrapolation across the saddle point was carried out for the velocity, 
density, and temperature. The individual species were also extrapolated, but 
the downstream calculations appeared to be sensitive to small changes 
("errors") in the species concentrations; hence, for the results shown, only 
the velocity, density, and temperature were extrapolated, and the species 
concentrations were chosen to be those in exact equilibrium with the 
extrapolated temperature and density. 

Sketches (k) and (&) also show the result of "bad" extrapolation; that 
is, the circled points are the results of calculations made with deliberate 
errors in the extrapolations for p, u, and T (these correspond to the 
example of the curve marked CC ! in sketch (a)). In this case, the conver- 
gence of the integral curves to the saddle curve is evident. No cases to the 
contrary were found, although no "extreme" errors were tried. 
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Sketch (k) . - Integral curves near the saddle 
curve for a real gas in nonequilibrium flow 
through a nozzle; reservoir conditions, 

T Q « 10,000° K, p Q » 4470' atmospheres . 



Sketch (m).- Integral curves near the saddle 
curve for a real gas in nonequilibrium flow 
through a nozzle; reservoir conditions, 

T 0 10,000° K, p 0 « Ml -70 atmospheres. 



Sketch ($,).- Integral curves near the saddle 
curve for a real gas in nonequilibrium, flow 
through a nozzle; reservoir conditions, 

T q = 10,000° K, p Q = 4470 atmospheres. 

The eigenvalues for the two 
lower pressure reservoirs, given in 
equations (28b) and (28c), as well as 
the solutions for the nonequilibrium 
flow through the nozzle region under 
these conditions, are shown in fig- 
ures 3 through 6. The parasitic 
eigenvalues tend to disappear as the 
reservoir pressure drops, and, in 
fact, in the lowest pressure case 
they have, for practical purposes, 
vanished. The latter case, therefore, 
can be integrated all the way from 
x — -1 using the explicit method 
1.6 defined in equation ( 20 ) . 

Notice that for the conditions 
given in equation 28(c) (low reser- 
voir pressure) the nozzle is out of 
equilibrium well upstream of the 
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minimum section. In this case, all of the species, as well as the velocity 
and thermodynamic variables, had to be extrapolated across the saddle point. 
For these conditions, however, this procedure appears to give quite satisfac- 
tory values as the results shown in figure 6 indicate by their smoothness. 


Calculations Made Near 1 Equilibrium 

As mentioned in the introduction to this section, some special problems 
arise when the flow is nearly in equilibrium, due to a numerical indeterminacy 


in the species production Q 1 (see eqs . ( 4 )). This 
the product of X^, the "degree of nonequilibrium," 


librium is approached, 
is given by 


For example, a typical Xj_, 


o' — 
- 1.0 


Sketch (n).- The fluctuation of du/ds in 
calculations made very near equilibrium; 
upstream from nozzle. 


Sketch (o).- The fluctuation of du/ds in 
calculations made very near equilibrium; 
downstream from nozzle. 


indeterminacy arises in the 
(ref. 7 ) Xi -* 0 as equi- 
N 2 2 2N, 

( 32 ) 


for the reaction 

1 - 


X P = 


P 7g* 




where K 2 is the equilibrium constant 
for the reaction specified. 

All computations were made in 
8 -place (27 -bit) floating-point arith- 
metic, and in the numerical calcula- 
tions under near equilibrium 
conditions the subtractions needed to 
form the X-j_ (e.g., eq. (32)) lost 
some of the significant figures. 
Although a quantitative analysis of 
the processes involved has not been 
carried out, the numerical evidence 
indicated that the loss was not disas- 
trous. It is thought that this is due 
to the fact that the product of step 
size and driving eigenvalues was suf- 
ficiently small that valid information 
was obtained from whatever significant 
figures remained. At any rate, suffi- 
cient information was retained to 
provide the accuracy indicated in 
figure 2 for the high-pressure reser- 
voir conditions (28a) and the flow 
model implicit in equations (5) • 

In near equilibrium calculations, 
most of the derivatives fluctuated 
within a band on successive steps . 

For example, typical results for 
du/ds on the upstream and downstream 
side of the throat are shown in 
sketches (n) and (o) . This kind of 
behavior for the derivative terms is 
typical of the unconditionally stable 
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implicit methods when they are being used in the presence of large negative 
eigenvalues. The following is offered as an explanation: 

1. The solution to the difference equations (in linearized form, at 
least) depends upon the sum of a number of terms represented (for the implicit 
modified Euler method) by 

— 1 ~~| n 

1 + - ha k 

° k — hr 

1 - g ha k_ 

where n is step number, cr k is an eigenvalue in [A n ], and C k is a constant 
that depends on the initial conditions. 

2 . If ha k is negative, the magnitude of the term in brackets is always 
less than 1, but if hcr k is a large negative number, the magnitude of the 
term within brackets is not much less than 1, and, furthermore, it is always 
negat ive . 

3 . In effect, then, the solution of the difference equation is composed 
of two kinds of terms. For example, 

driving terms 

f a — — — ^ 

w = c 0 (o .021 . . .) n + qf -0.035 . . .) n + . . . 


parasitic terms 

— — ^ ^ 

+ C k ( - 0.999 • • -) n + c k+1 ( -0.999 . . -) n + . . . 

b. The terms (- 0-999 ■ • -) n oscillate between +1 and -1 as n is 
even or odd and if the coefficients C k , C k+1 , • • • are not exactly zero 
(they cannot be exactly zero in truncated floating-point arithmetic) . This 
will appear as "noise” at some level of significance in the numerical 
solution. 

5. Apparently, for the mathematical model chosen and for the arithmetic 
used to study it, this noise is just below the level of usable information in 
the values of the dependent variables at near equilibrium conditions. 

When calculations were made even slightly away from equilibrium, the 
fluctuations (to the scale shown) disappeared and the numerical calculations 
appeared "normal.” For example, the transition from fluctuating to relatively 
smooth values of the derivative du/ds is seen in sketch (o) to occur about 
1 cm after passing the minimum section in the high pressure case. To the 
scale used in figure 1, the gas is not yet out of equilibrium until about 3 cm 
downstream from the throat. 
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CONCLUDING REMARKS 


For problems that require the numerical integration of differential 
equations through a saddle-point type of singular ity, it is necessary (l) to 
identify the location of the saddle point and then (2) to use special numeri- 
cal techniques to obtain the solution of the differential equations in the 
saddle -point vicinity. If the eigenvalues of the system having different 
sign are also of different magnitude, the technique of backward integration, 
that is, integration from the saddle-point vicinity outward, may be useful. 

Large parasitic eigenvalues may be present in the numerical integration 
of the differential equations governing the flow of a gas out of chemical 
equilibrium. In the model chosen, this was particularly true near the throat 
of the nozzle, and the severity of the parasitic behavior increased with 
increasing density. Use of an implicit numerical method is recommended for 
the most severe cases. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035 , Jan. 6, 1969 
129-01-02-05-00-21 
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(a) Thermodynamic variables* 

Figure 2.- Solution for the flow through a nozzle of a real gas in 
nonequilibrium; T 0 = 10*000° K, p = ^70 atmospheres. 
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Figure 2.- Concluded. 
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Figure 3*- Local eigenvalues in the equations representing the flow through a nozzle of a real gas in 
nonequilibrium; reservoir conditions, T = 8000° K, p - 100 atmospheres. 
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(a) Thermodynamic variables. 

Figure 4.- Solution for the flow through a nozzle of a real gas in nonequilibrium; T 0 = 8000° K, 

p = 100 atmospheres. 
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(a) Thermodynamic variables. 

Figure 6.- Solution for the flow through a nozzle of a real gas in nonequilibrium; T 0 = 10,000° K 

p =0.2 atmosphere. 
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Low density case 
T 0 = 10,000° K 
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(b) Species concentrations . 

Figure 6.- Concluded. 



NASA -Langley, 1969 19 A -310 4 



National Aeronautics and Space Administration 
Washington, D. C. 20546 


POSTAGE and pees paid 
NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 


OFFICIAL BUSINESS 


FIRST CLASS MAIL 




POSTMASTER: 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 


"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the tv ides t practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and Notes, 
and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 


